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Abstract: This paper compares three boundary conditions, i.e. transmitting boundary condition, Sarma 
absorbing boundary condition and the uniaxial complete matched layer absorbing boundary condition for 
simulation of ground penetrating radar (GPR) by the time domain finite element (FEM) method. The 
formulations of the three boundary conditions for the FEM method are described. Their effectiveness in 
absorbing the incident electromagnetic waves are evaluated by the reflection coefficient on the boundary 
of a simple GPR model. The results demonstrate that UPML boundary condition can yield a reflection 
coefficient smaller than -50 dB, which is -20 dB smaller than other two boundary conditions. 
Keywords: Ground penetrating radar (GPR); uniaxial complete matched layer absorbing boundary 
condition; Sarma boundary condition, transmitting boundary condition; finite element method 

I. Introduction 

Ground penetrating radar (GPR) has been widely applied in many fields. Numerical simulation plays a 
great role in the interpretation of GPR data and the inversion of physical parameters. Finite element method 
(FEM) is widely used due to its easy adaption to a complex subsurface dielectric model with a triangle mesh. 
Shen et al [9] have performed a 2D forward GPR modeling by substituting the electromagnetic wave 
transmitting equation with the sound wave equation; Di and Wang [10] introduces the complex dielectric 
permittivity into the finite element equation using the Galerkin method. Tiao et al [ll]perform a GPR numerical 
simulation on dispersive media using a discrete time-domain method. Arias et al [12] applies the GPR finite 
element technique to the imaging of archaeological remains. For GPR numerical simulation, absorbing 
boundary condition is indispensable for eliminating the undesired boundary reflection. An effective absorbing 
boundary condition is critical for improvement of the accuracy of GPR numerical simulation. Many research 
efforts have been paid to the absorbing boundary conditions for the GPR forward modelling using the finite 
different time domain (FDTD) algorithm. There are several absorbing boundary conditions working with 
different Attenuation mechanisms. For an example, the radiant boundary condition (Bayliss and Turkel E. 2009) 
and absorbing boundary condition are based on the separate wave equation(Engquist and Majda. 1977); ultra- 
absorbing boundary condition (Mei. 1992), and complete matched layer(Gedney et al. 2001; Berenger. 1994,) 
absorbing boundary condition have been applied. For the GPR simulation by the FEM method, Di and Wang 
(2005) use a transmitting absorbing boundary condition; Sarma absorbing condition has been used in GPR 2D 
numerical simulation by Feng (2013). UPML absorbing boundary condition also has applied in GPR 2D 
dispersive medium numerical simulation (Tiao L, 2006). In this paper, the principle and the formulations of the 
transmission absorbing, Sarma absorbing and UPML absorbing boundary conditions are going to be presented. 
Their performances while applied to the FEM simulation are tested. 

II. Fern Applied To Wave Equation Of Gap 

The electric field radiated by a GPR source is governed by wave equation, as expressed by (Di, et al., 1999) 

dt /us s y L > 

where £ is the dielectric permittivity, /J is the magnetic permeability, <J is the electrical conductivity, E is 
electric-field intensity and S E is the electric field source. We divide the 2D region into rectangular elements and 
let E be a linear function of position coordinates in each element, as given by 

E = N T E e (2) 

where is the shape function for linear interpolation and E e is the column vector whose components are E . 
at the node / of the element where i = 1, 2, 3, 4 . Using the Galerkin method, the general 2D finite-element 
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equation in time and spatial domain can be expressed as 

M^K&+KE = S E (3) 

where M is the mass matrix, K' is the damping matrix, K is the stiffness matrix, & anc [ ffi are the first- and 
second-order derivatives of electric field, respectively .We can adopt the central difference method in time 
domain to solve the Eq. (3). 

III. Absorbing Boundary Conditions 

In this section, we apply the transmitting boundary condition, sarmaboundary condition and UPML 
absorbing boundary condition to the FEM simulation of GPR data. 
3.1 Transmitting boundary condition 

The transmitting boundary condition is established by directly simulating the common kinematic properties 
of various one-way waves. The principle is that the one-way waves are transmitted through the boundary 
interface at one point on the artificial boundary. This means that the one-way waves can be expressed as the 
superposition of a series of outgoing plane waves, because there is no specific limit to the speed at which these 
plane waves transmit along the artificial boundary, nor any limitation on their shapes. Assuming that all the one- 
way waves are transmitted along artificialboundary at the same velocity, we obtain a formula for the 
transmitting boundary condition [24-25]. 
Suppose 

d 2 E 1 2 cr dE 

LE = — V 2 £ + = V(4) 

dt jus s dt 

At each element, the electrical field E can be expressed as 

E = ^N { E. ' ( 5 ) 

i=i 

where E i and N. refer to the electric field and the shape function of the node I in each element. The residual 
quantity R can be expressed as 

R — LE—S E . (6) 

According to the Galerkinfinite element method, the 2D finite element equation can be expressed as 

N T rdxdy = 0, (7) 

e 

where Q refers to the area of the element and Y is the vector expression of R According to the paraxial 

approximation deduced by Clear.bout[26], the down-going, up-going, left-going, right-going wave equations in 

a homogeneous isotropic medium, arerespectively expressed as 

dE 1 dE dE 1 dE dE 1 dE dE 1 dE 

— = , — = , — = , — = , (8) 

dt v dy dt v dy dt v dx dt v dx 

According to the definition of the normal derivative, we know that 

d d d 

— = — n x + — n , (9) 
dn dx dy 

where ^ x and n y refer to the direction cosines of the boundary external normal.By importing Eq. (5) and (6) 
into Eq. (7) and applying Gauss theory, we can obtained 

^ r d ( dE T ^ r dE T r dE T r dE T r dE T 

Y\ — — n dxdy=\ — Ndr+\ — Ndr+\ — Ndr\ — Ndr, (io) 

p v Jr, Q n Jr d Q n )r r Q n Jr M Q n 



dx\ dx 



J 



where 1^ , Y r , Y d , Y u representthe left, right, bottom, upper boundary respectively. Importing the boundary 
conditions of Eq. (9) into Eq. (10) we can obtained 
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e dE T _ f 8E T _ f DE T _ r DE T 
NdF+\ NdT+\ NdF+\ 

Jr ' dn Jr d dn Jr ^ dn Jr » dn 



NdF 



= v$x 



(11) 



l t N,N 2 (n x + n y ) dY + N 2 N 3 (n x + n y ) dT 

+l r N i N * (n x +n y )dT + N,N 4 (n x +n y )dT 

where the first-order derivative of the electric field in time, and v refers to the transmission velocity of 
electromagnetic waves of medium. Therefore, the damping matrix of the boundary F can be expressed as 

" =!.". =1.(1 d -f NdT+ l d ~f NdT+ l d 4- NdT+ l d -ffdT )m 

q q 1 Sn JT d on :r r an Jr « on 

By importing boundary damping matrix into the GPR finite element equation, we can obtained 

Ml l ^(K , + F)$+KE = S E . (13) 

3.2 Sarma absorbing boundary 

The Sarma absorbing boundary is a kind of attenuating absorbing layer. The basic principle is that an 
absorbing layer is added the boundaries of the simulation model to absorb the incident GPR waves. The medium 
of the absorbing layer is a kind of high attenuation characteristic media. The attenuation coefficient of the 
absorbing layer should be built on reasonable physical theory, and expressed as function of media parameters. 
Sarma (Sarma et al., 1998) proposed an absorbing boundary condition by adding an attenuating layer of a 
certain thickness around the boundaries area according to the material property to estimate the constant 
proportion coefficient. Figure 1(a) shows a sketch of the Sarma boundary condition. Rayleigh gives a classic 
method to calculate the damping matrix F .This method uses the whole mass matrix M and total stiffness 
matrix K to calculate damping matrix F .The matrix F can be expressed as 

[F] = a 0 [M] + a\K],{U) 

where a Q , a x can be expressed as 

_ 2co i co ] {^co ] _ 2(^0). -Zq) 



,(15) 



CO - CO 



CO - CO 



The quantities, CO. and CO- are the inherent radiosfrequency of the / and j media, and and are the 

corresponding damping ratios, respectively. We apply the Sarma absorbing boundary to process the reflected 
wave. We substitute the matrix F into equation (3).Then EFM equation in the absorbing boundary area can be 
obtained as 

M^(K' + F)&+KE = S, (16) 
In the absorbing layers areas, the proportionality coefficient a 0 and a x can be calculated according to equation 
(15). 

3.3 UPML absorbing boundary condition 

According to the electromagnetic theory, Maxwell equations in a UPML medium can be expressed as 

VxE = -}cDjuS-H (17) 

WxH = ]cosS'E{\%) 

where S is a parameter which expresses the characteristic of an anisotropic medium. It can be expressed as 
follows in two dimensions. 

0 0 



S = 



where s t = 1 + cr- /( j cos 0 ), i = x, y 

To solve rotation on both sides of Eq. (17), and put Eq. (18) into it, we can obtain 



s y l s x 
0 
0 



s x /s y 



0 



(19) 



-S~ l VxVxE = co 2 sS-E (20) 
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Substitute Eq. (19) into Eq. (20) and we obtain 



{\ + a x li]cos)Y 
0 
0 



0 0 

(l + a y /(]co£)) 2 0 
1 

\2 



oYs- 



0 

'<\ + a y l<$a>e)Y 



VxVxE= 



0 



0 
0 



{\ + cj x I{]cos)Y 
0 



0 
0 

(\ + cj x IQcos)Y(\ + cj J(]cos)Y 



(21) 



E 



Apply Fourier transform to the above equation, and we can obtain 



1 ^ „ 2 T dVxVxE 1 x2 dVxVxE 
— VxVxE + / + — -I 

.2 ^2 



//^ 0 //6*q dt z 

Applying theGalerkin method [15] to Eq. (22), the FEM equation in a UPML medium can be expressed as 

M,#f M 2 £ + M 3 J' Edt + M 4 J' Edtdt -KE- K x &- K 2 $= S (23) 
The express of the matrixes in Eq. (23) is as follows 

M = f N T Ndn 

Je 

K=±j 

iip Je 



JUS ' 





'aw' 


^ + dN 


(8N^ 




dx 1 

V 




I + &y 




) 



dQ 



" 'I 



flS Q £ 



dNi 


'dN~ 


f dN 


fdN^ 




dx 1 
V 


v 3x , 


I + dy 




J 



dQ 



M 2 =-^- I {J 2 +K 2 +2-H)\ e N T NdQ 



M £ 0 2 s 3e 



'dNi 


'aw' 


^ + dN 


(dN^ 




dx \ 

V 




) + dy 




) 



dQ 



M 3 = — 



N T NdQ 



K H\ 

j ( 

M 4 =-^-Hj N T NdQ 
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IV. Comparison test and results 

To illustrate absorption effectson the incident waves at the artificially truncated boundary using different 
absorbing boundary condition, a lOmx 10m homogeneous model was established, with a pulse excitation 
source located at the center of the simulation area, as shown in Fig. 1. The wavelte of the GPR pulse excitation 

source is f(t) — t 2 e~ at sin OJ^t , where a> 0 is the central frequency of the transmitting antenna (in this example, 

#; 0 =100 MHz). The attenuation velocity is governed by the coefficient a , which in this case equals to 0.93co 0 - 

The sampling interval time is 0.1 ns. The dielectric constant of the medium is G — 6.5 , the conductivity 
<J = 0.002 S/m, and the mesh element size = 0.0025 m 2 . Absorption effects resulting from different absorbing 
boundary conditions were observed by wave-field snapshots at various times. Fig. 2. Fig. 3 and Fig. 4 show the 
snapshots of wave-field surrounded by transmitting boundary condition, Sarma boundary condition and UPML 
boundary condition respectively. In Fig. 2, we can see that the GPR wave reaches the edge of the area at 32 ns, 
and the reflected wave begins to form at the artificial truncated boundary. In the snapshot at 38 ns, we can see 
the GPR wave energy reflected strongly by the artificial truncated boundary, greatly affecting the target area. 
This also indicates the necessity for effective processing of artificial truncated boundaries. Fig. 3 and Fig. 4 
show wave field snapshots at 12, 32 and 38 ns by using the transmitting boundary condition and Sarmaboundary 
condition respectively. As seen in the figures, wave energy again reaches the boundary area at 32 ns. By 38 ns, 
however, the majority of the electromagnetic wave energy has been transmitted through the boundary, with only 
a small portion reflected back. By comparison to the energy of the strongly reflected wave on the truncated 
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boundary at the same time (Fig. 5(c)), we see that the transmitting boundary condition and Sarmaboundary 
condition has a remarkable effect. However, wave field snapshot of 32 ns shows that there is still some reflected 
energy from the truncated boundary. Compared with Fig. 3 and Fig. 4, wave filed snapshots of 32 ns in Fig. 5 
which is by using UPML absorbing boundary condition have perfectly absorbing effect, the reflected wave are 
nearly absorbing. 

Fig. 6. shows the recorded signal and their envelope at point A in Fig. 1 under different absorbing boundary 
conditions. The recorded signal contains the transmission signal and the reflected signal from the boundary. 
From these figures we can see that under the UPML boundary condition, the reflected energy is least. The 
refection efficient of no absorbing boundary condition , transmitting boundary condition, sarmaboundary 
condition and UPML boundary condition calculated after a spreading correction ofthe wave front are -9.68 dB, - 
24.8 dB, -34.7 dB,-55.8 dB respectively. Compared with other two absorbing boundary conditions, the UPML 
boundary condition has best absorbing effect. 
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Rx 
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A 


B 



absorbing layers 



Fig.l. Sketch drawing of the Sarma boundary condition model 




(a) 12 ns snapshot(b) 32 ns snapshot (c)38 ns snapshot 
Fig. 2. Snapshots of wave field without boundary conditions 




(a) 12 ns snapshot (b) 32 ns snapshot(c) 38 ns snapshot 
Fig. 3. Snapshots of wave field with the transmitting boundary condi 
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(a)32 ns snapshot (b) 38 ns snapshot(c) 46 ns snapshot 
Fig.4. Snapshots of wave field with Sarma absorbing boundary condition 



f 5.00 1 



o 



^^^^ 

L™ ^ 1 4D0^ J" 



1.25 2.50 3.75 5.00 6.25 7.50 8.75 10.0C 



1.25 2.50 3.75 5.00 6.25 7.50 8.75 10.00 



1.25 2 50 3.75 5.00 6.25 7.50 8.75 10.00 



(a) 12 ns snapshot (b) 32 ns snapshot(c) 38 ns snapshot 
Fig. 5. Snapshots of wave field ts with the UPML boundary condition 
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Fig. 6. Incident and reflected signal from different absorbing boundary, (a) No absorbing boundary condition;(b) 
Transmitting boundary condition; (c) Sarma absorbing boundary condition and (d) UPML absorbing boundary 
condition. 
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IV. Conclusion 

Transmitting boundary condition, Sarma absorbing boundary condition and UPML absorbing boundary 
condition for simulation of ground penetrating radar (GPR) by the time domain finite element (FEM) method 
are described and compared by the reflection coefficient.The results demonstrated that UPML boundary 
condition can yield a reflection coefficient smaller than -50 dB, which -20 dB smaller than other two boundary 
conditions. 
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